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SUMMARY 

We develop a Discrete Element Method (DEM) for elastodynamics using polyhedral elements. We show that for 
a given choice of forces and torques, we recover the equations of linear elastodynamics in small deformations. 
Furthermore, the torques and forces derive from a potential energy, and thus the global equation is an Hamiltonian 
dynamics. The use of an explicit symplectic time integration scheme allows us to recover conservation of energy, 
and thus stability over long time simulations. These theoretical results are illustrated by numerical simulations of 
test cases involving large displacements. Copyright © 2010 John Wiley & Sons, Ltd. 
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1. Introduction 

Particle methods are meshless simulation techniques in which a continuum medium is approximated 
through the dynamics of a set of interacting particles. Two main classes of particle methods can be 
distinguished : Discrete Element methods (DEM), which rely on the contact interaction of material 
particles by means of forces and torques, and Smooth Particle Hydrodynamics (SPH) methods, in 
which the continuum is discretized by localized kernel functions. 

Discrete Element methods consist in the resolution of the equations of motion of a set of particles 
submitted to forces and torques. It is thus possible to account for a variety of phenomena (behaviour 
laws, models, scales,...) using a single numerical method. A wide variety of Discrete Elment methods 
have been designed changing the expression of the forces, with particular attention devoted to specific 
aspects. Discrete Element methods have first been developed by Hoover, Arhurst and Olness J20) in 
models for crystalline materials. Their application to geotechnical problems was carried out by Cundall 
and Strack [4|, and their use in granular materials and rock simulation is still widespread J36][37) . 
Discrete Element Methods have also been used to simulate thermal conduction in granular assemblies 
p0| or fluid-structure interaction p6| . The model is also able to account for grain size effects ]21) , and 
to treat fracture in a natural way. Discrete Element methods used for granular materials generally 
describe particles as spherical elements interacting via noncohesive, frictional contact forces [37) . 
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For brittle materials, models also use unilateral contact forces, combined with bonds which simulate 
cohesion |36|. Kun and Herrmann developed a combination of the contact model with a lattice model 
of beams to account for the cohesion |26|, which has been extended to Reissner models of beams 
to simulate large rotations of the material [5 21]. The authors use Voronoi tesselations to generate 
the polygonal particles. However, the results obtained still depend on the size of the discretization 
(which physically corresponds to the size of heterogeneities) pT) . The effective macroscopic Young 
modulus and Poisson ratio highly depend on the isotropy of the distribution of the particles and are 
only empirically linked to their microscopic value for the Reissner beams p6) . 

In a different approach, SPH methods describe the particles as smooth density kernel functions. The 
kernel functions are an approximation of the partition of unity. The continuous equations of evolution of 
the fluid or solid material therefore induce the dynamics of the particles. Originating from astrophysical 
compressible fluid simulations fT2"}[33) , SPH was extended to incompressible fluids [35 1 and to elastic 
and plastic dynamics p2| , and used for fluid-structure interaction with both domains discretized with 
SPH [2|. A state of the art review of the method with applications to solid mechanics is presented 
in p9||7 SPH preserves the total mass of the system exactly. However, in tensile regime, unphysical 
clusters of particles tend to appear in situations where a homogeneous response is expected |40|. 
Hicks, Swegle and Attaway advocate the smoothing of the variables between neighbouring particles to 
stabilize the method, rather than introducing artificial viscosities 1 18 1. Bonet and Lok have addressed 
the issue of angular momentum preservation, and show that rotational invariance is equivalent to 
the exact evaluation of the gradients of linear velocity fields, which can be achieved either through 
correction of the kernel function or through a modification of its gradient pj. In order to circumvent 
the difficulties affecting SPH, Yserentant developed the Finite Mass method, in which particles of 
fixed size and shape also possess a rotational degree of freedom (spin). The method achieves effective 
partition of unity, and thus preserves momentum, angular momentum and energy, ensuring stability 
g2). 

The Moving Particle Semi-implicit (MPS) method is a variant of the SPH method developed by 
Koshizuka. It consists in the derivation of the dynamics of a set of points from a discrete Hamiltonian 
J23) . As in the SPH method, the differential operators are approximated by a kernel function of 
compact support. The expression of the approximated differential operators is inserted in the classical 
Hamiltonian of the system, and by application of Hamilton's equations, the dynamics of the discretized 
system is obtained. To preserve the Hamiltonian structure of the dynamic of the system through time 
discretization, the authors use symplectic schemes (39). The MPS method has been used initially for 
free-surface flows [23,24], and has been extended to nonlinear elastodynamics [25 39 1 and to fluid- 
structure interaction [29 1. Using similar ideas, by deriving the dynamics of the system from a discrete 
Hamiltonian, Fahrenthold has simulated compressible flows f22") and impact events with breaking of 
the target (8]|9j. 

These methods show the importance of the preservation of momentum and energy for the accuracy 
and stability of the scheme over long-time simulation. The use of symplectic schemes ensures the 
preservation of the structure of Hamilton's equations by the numerical time integration, and therefore 
the preservation of momentum and energy [15]. Simo, Tarnow and Wong note, however, that while 
ensuring the stability of the simulation for small time steps, the symplectic schemes fail to preserve 
exactly energy and become unstable for larger time steps [38 1. They derive a general class of implicit 
time-stepping algorithms which exactly enforce the conservation of momentum, angular momentum 
and energy. The algorithms are built in order to preserve linear and angular momentum, and energy 
conservation is enforced either with a projection method (projection on the manifold of constant 
energy) or with a collocation method. The algorithm is used for nonlinear elasticity in large deformation 
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using finite element methods 1 13 28 , 38 1 and for low-velocity impact [ 17]. 

In this article, we extend and analyze the Discrete Element method initially introduced by Mariotti 
[34| . Combining a Discrete Element Method with a lattice model of beams, we are able to account 
for the cohesion of the material, and analytically recover the macroscopic behaviour of the continuous 
material. The method, Mka3D, has been successfully used to simulate the propagation of seismic waves 
in linear elastic medium p4| . Here, we extend the properties of the algorithm to the case of large 
displacements without fracture. Contrary to usual Discrete Element methods, we are able to derive the 
microscale forces and torques analytically from the macroscopic Young modulus and Poisson ratio, and 
to prove the convergence of the method as the grid is refined. In addition, as in MPS methods, we derive 
the forces and torques between particles from a Hamiltonian formulation. Using a symplectic scheme, 
we ensure the preservation of energy over long-time simulations, and thus stability of the method. 
This allows for the simulation of three-dimensional wave propagation as well as shell or multibody 
dynamics. The paper is organized as follows. In section [2] we describe the lattice model used. We 
introduce the Hamiltonian of the system and we derive the expression of forces and torques chosen to 
simulate linear elasticity. In section[3] we show that these expressions lead to a macroscopic behaviour 
of the material equivalent to a Cosserat continuum, with a characteristic length of the order of the size 
of the particles. Hence, the model is consistent with a Cauchy continuum medium up to second-order 
accuracy, in the case of small displacement and small deformation. The microscopic values of Young 
modulus and Poisson ratio yield directly the macroscopic values, and we can choose Poisson ratio in 
the whole interval (—1, 0.5). In section|4] we then describe the symplectic RATTLE time-scheme ]T5) , 
which allows us to preserve a discrete energy over long-time simulations. These theoretical results are 
illustrated by numerical simulations of test cases involving large displacements in section [5J 



2. Description of the method 
2.1. Geometrical description of the system 

In order to discretize the continuum material, several methods have been suggested for Discrete 
Element Methods. Most authors working on granular materials use hard spheres, in order to simplify 
the computation of contacts between particles, as the exact form of the particles is mainly unknown. 
However, in the case of the simulation of a continuous material, this method is not adapted as the 
interstitial vacuum between spheres is inconsistent with the compactness of the solid. In addition, the 
difficulty to obtain a dense packing of hard spheres, and the problem of the expression of cohesion 
between the particles, have led us to use Voronoi tesselations instead, as suggested in 15] [26). The 
particles are therefore convex polyhedra which define a partition of the entire domain. As we shall see, 
this method allows us to handle any Poisson ratio v strictly between —1 and 0.5, independently from 
the size of the particles. On the contrary, most granular sphere packing methodologies account for a 
limited range of v, which is size dependent. 

The following parameters are relevant to describe the motion of a given particle I : X_j and v_j denote 
respectively the position and velocity of its center of mass (v_j = < ^ L ), Q denotes the orthogonal 
rotation matrix of the frame attached to the rigid particle, and the angular velocity vector f2 7 is uniquely 
defined by : 

dQ 

m = -=% r , (i) 
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Particle I 



+Xj 



Particle J 



Figure 1 . Geometric description of the particles 

where the map j : M. 3 — >• R 3x3 is such that : 

Viel 3 , Vy G M 3 , i(x) -y = xAy 

Finally, the material of particle / is described by its mass mj, its volume Vj and its principal moments 
of inertia I}, ij and if. We suppose the local frame attached to the particle is attached to the principal 
axes of inertia (e}, ef , e 3 ). The matrix of inertia in the fixed frame is given by : 



R 



(2) 



where Rz is the matrix of inertia RP written in the inertial frame 



1? 



We also define the parameters dj, dj and d 3 - as : 

d 



1} 






















If 



1} + Ij + 1] 



Pi, 



1,2,3 



and we introduce the following matrix D defined in the inertial frame : 

/ d) 
£ = 
\ d 3 j 

The Discrete Element Method relies on the computation of forces and torques between nearest 
neighbours particles. We denote by Vi the list of the neighbouring particles linked to particle /. For 
each link between two particles / and J, we define Pjj the center of mass of the interface, Su the 
surface of the interface, the distance between particles I and J : 



D u = \\XjXj 
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and the initial exterior normal vector for link IJ : 

We define two normalized orthogonal vectors of the interface s ZJ and i J7 = n J7 A s 7J , serving as 
references to evaluate the torsion between particles / and J. 

These parameters are given a fixed value at the beginning of the computation. D\j and rttjj 
respectively denote the initial values for Djj and njj. The particles are therefore assumed to be 
rigid. However, compressibility effects are taken into account through the expression of interaction 
potentials. 

In addition, we define the following quantities : 

• the displacement at the interface between particles / and J : 

Auu = Xj-X l + Q j - XjPu - Q x ■ X^Pu 

• When particle / has several free interfaces (i.e. not linked to another particle), these surfaces are 
marked as stress-free. To account for the free deformation of the particle in these directions, free- 
volume V} is defined as the sum of the volumes of all pyramidal polyhedra with a free surface 
as basis and X® as summit. 

• the volumetric deformation e\ of particle / is defined as the sum of all contributions of the 
deformations of the material links of particle /. We have assumed that the bending of the link 
between two particles does not affect volume, as long as the centers of the interface of the two 
particles stay in contact. The corrective term on the volume is active only on particles having a 
free surface, and accounts for the boundary condition q_ ■ n = 0. We derive it in Appendix |I] 

• The interpolated volumetric deformation for link (IJ) : 

eh = \{e v i+e v j) 
2.2. Expression of the Hamiltonian of the system 

We denote by E the Young's modulus and by v the Poisson's ratio for the material. The Hamiltonian 
formulation of the elastodynamic equations on a domain £1 is as follows : 

H («>P)= f ^P'P+ U (£) < 3 ) 

where q is the displacement field and p = pv is the density of momentum. U(q) is the potential 
energy of the system. It can be expressed in terms of the stress tensor a_ and the linearized strain tensor 

i = |(y9+y? T ): 

U(q) = W(g) = i / g(g) : £ (4) 
In the case of Cauchy linear elasticity, we use the constitutive relation 



E Ev 
^ ^ —u £ = + (1 + „)(!- 2.) tr ^ (5) 



Copyright © 2010 John Wiley & Sons, Ltd. 
Prepared using nmeauth.cls 



Int. J. Numer. Meth. Engng 2010; 0:0-0 



AN ENERGY-PRESERVING D.E.M. FOR ELASTODYNAMICS 



5 



to derive the expressions of W(e) and U (q) : 

g w - a /„ = a + 2(1 + ^ _ a) 3 ro 

We choose to discretize the Hamiltonian formulation as a discrete Hamiltonian if/j. The 
displacement field q is derived from the values of (Xj, Q ). The density of momentum derives from : 

Tj = mivx (8) 
g I =Mi)-g l -B zl (9) 



We define : 



H h (X, g,T,Q = \j2 ~Tj ■T I + l -Y j ir(g I -R- 1 ■ £/) + U h (X, £) (10) 

The discretized potential energy is split into three terms : 

u h {x,g) = U t (X,g) + U d {X,g) + U f (Q) 

Ut(X, Q) corresponds to the first term of (|6ji : we approach the strain of the link (IJ) in the direction 
Hu 1 ' Rij by the normalized displacement jjn-Aujj, and we use the approximation : 

|:£« E (i-Zuf (ID 



We therefore write 

u t (x,g) - ^ Sjj— — ^- — 

This energy accounts for the deformation of each link between two particles. 

Ud(X_, Q) corresponds to the second term of doll : we approach the trace of the strain tr(e) in particle 



I by the sum of the normalized displacements ef for links surrounding /. A corrective term is added 
for cells having a free boundary : 

This energy accounts for the global volumetric deformation of each particle. 

The former two terms are sufficient to recover the equations of elastodynamics inside the solid. 
However, for the method to be able to cope with thin one-element shells, we add the pure flexion term 

Uf(Q) : 



u f(Q) = - E §r (««(£, • ak) ■ (£ z ■ oSj) +«,(£ J • sjj) ■ (q i ■ sjj) + a t (q j - 1„) • ■ t, 
(u) IJ 



This term accounts for the flexion between particles. The coefficients a n , a s and at are chosen to 
recover the exact flexion and torsion of a beam, and are detailed in Appendix [IT] 
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2.3. Derivation of the forces and torques between particles 



We use Hamilton's equations for the system ( 10 1 : 



Ki = 5^ d2) 



dTj 



(13) 



£' = -!; 

^ = -ag; + ^-£, (15) 

where A 7 is the symmetric matrix of the Lagrange multipliers associated with the constraint Q T -Q = 

M- ' ~ ! _/ 

Equations ( p"2} and (\3\ give us the usual kinematic relations between position and velocity : 

Xj = mj 1 ^ = Vj 

The derivation of forces and torques from the potential energies is carried out in Appendix [ill] We 
obtain mjij = F_j j where F_u, the force exerted by particle / on particle J, is given by : 

Su E „ Ev „ ( 1 . 1 

-\mku-ikum 

(16) 



ElJ = Dfj " + SlJ {l + v){l-^f IJ + ^J— IJ " ^7 ( — " ■ - Ij) - U 



This expression can be seen as a discrete version of Hooke's law of linear elasticity 

E Ev 



—tr(e)Id (17) 
= l + i/= (l + v)(l-2v) y = J = 

using the previous analogies between jyn-Aujj and e, e| j and tr e, and noting that a ■ n is a force per 
surface unit (a pressure). 

For the rotational part, we define the two following torques : 

— u = w J ^M/^ A ^» + (i + ,Ki- 2 ,) g ^ ( a • ^ A ^ (18) 

Ml j = {a n {Q i ■ n° u ) A {Q j ■ n°j) +a s (g i ■ s LJ ) A {g j ■ s u ) +a t (Q j • t LJ ) A (Q j ■ t u j) 

(19) 

We note the fact that M^j corresponds to the torque at the center of mass of the force Fjj exerted by 
particle J on particle / at point Pjj : 

M t u = (Q I - X° I P IJ )AF IJ 
and Mjj is the flexion-torsion torque. We get the equation on the angular velocity : 

JeVi 

In the case when exterior forces and torques are applied to the system, they are to be added to the 
internal forces and torques computed above. 
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3. Consistency and accuracy of the scheme 



In this section, we investigate the consistency and the accuracy of the scheme. We first propose a 
modified equation for small displacements and small deformations. As the equations obtained are 
coupled dynamics for displacement and rotation, we compare the model with Cosserat generalized 
continuum, and recover a Cauchy continuum as the spatial discretization h tends to zero. 



3.1. Modified equation for the scheme 



The modified equation approach is a standard scheme analysis where a set of continuous equations 
verified by the approximate solution is seeked for. These modified equations should be an approximate 
version of continuous equations derived from physics. 

In order to be able to carry out a Taylor developments of the displacement, we place the points of the 
Voronoi tesselation on a Cartesian grid. The Discrete Element method can be seen, in this simplified 
case, as a Finite Difference scheme. 

We assume that no exterior force and no exterior torque are applied on the system. The displacement 
£ of particle / is given by : 

We assume that £ is a regular function on the domain, and we can therefore expand £ at point / 
with Taylor series if J £ V/. We denote Ax, Ay and Az the grid steps in each direction, and h their 
maximum. 

We assume displacements and rotations to be small. We denote 6 X , Oy and B\ the small rotation 
angles around axes x, y and z. 



Using ( 16 1, a simple Taylor development of the equations of motion yields for the displacement : 



pCx = 



+ 



E 



1 

Ev 



d% | d 2 £ x , 8 2 ij x 
dx 2 dy 2 
E fAx 2 d 4 ^ x 



12 dx 4 
Ax 2 d% 



dz 2 dy 

. Ay 2 d% 
12 dy 4 
Ax 2 d 4 £ y 



dO y \ | Ev f d 2 £, x | 8% | d% 

dz J (1 + v)(l — 2v) \ dx 2 dxdy dxdz 
Az 2 d 4 ^ x Ay 2 d 3 9 z Az 2 d 3 6 y 



(l + p)(l-2v) \ 3 da; 4 6 dx 3 dy 



12 dz 4 6 

Ax 2 d A £z 
6 dx 3 dz 



dy 3 6 
Ay 2 8% 
6 dxdy 3 



dz 3 



6 dxdz 3 



+ 0(h 3 ) (21) 



The same results hold for and £ z permuting the indices x, y and z circularly. 



Using ( 18 1 and ( 19 1, ((20b gives the equivalent equation for the rotation : 



Ay 2 + Az 2 
12 



E 



= _E_ fd^ _ diy 
1 + v \ dy dz 
Az 4 d b j y Ay 2 d 2 9 x 
120 dz 5 
Ay 2 + t 



4 dy 2 

d 2 e x 



12(l + i/) \dx 



26 x . 

6 

Az 2 d 2 e x 



Ay 2 d 3 i z Az 2 d 3 Cy Ay 4 d 5 i z 



12 dx 4 



dy 3 6 
Ay 4 d 4 6 x 



dz 3 120 dy 5 
Az 4 d 4 9 x 



8 dy 4 
- 2 'd 2 6 T 



Az 

+ T2~ \dy 2 
Ay 2 fd 2 9 x Az- 
U~ (lh 2+ T2 



:8 dz 4 
Ay 2 d 4 9 x 
12 dy 4 - 



dz 4 



0{h 5 ) (22) 
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The same results hold for y and 9 Z permuting the indices x, y and z circularly. 

We see that these sets of equations couple £ and 9, and by construction of the method, no constitutive 
law exists between £ and 8. The fact that a rotation remains in the equations can be compared to 
Cosserat continuum theory. We investigate this comparison in the following subsection. 

3.2. Comparison with Cosserat and Cauchy continuum theories 

In a Cosserat model for continuum media, the kinematics is described by a displacement field u and a 
rotation field <f>. A modified strain tensor e and a new curvature strain tensor k are introduced |7j : 

i = ym + in) 

We define t and fi the stress and couple stress tensors. We assume the following constitutive relations : 

t = Atr + fig + fi c g T (23) 

p, = crtr (re) Id + jk + /3k t (24) 

where A, fi, fi c , a, (3 and 7 are elastic moduli. 
The dynamical equations for the system are : 

pu = div t 

/ <f> = div p + e : t 

where p denotes the density, I is a characteristic inertia matrix, : denotes the double contraction 
product of tensors, and e is defined as follows : 

1 if (ijk) is an even permutation 
(&)ijk = i — 1 if (ijk) is an odd permutation 
otherwise 



Using the constitutive relations ( |23] > and p4| , the following equations can be obtained : 

pu = (A + p c ) Vdiv u + fiAu + (fi — /i c )c_url tf> (25) 
l c 4> = (a + /3)Vdiv <j> + 7A0 - 2(fi - p c )4> + (p - /i c )curl u (26) 
Identifying the terms of ( |25] > with equation pT) , we find : 

Si/ 

A 



(1 + i/)(1-2i/) 



l + i/ 
= 

and we therefore recover the classical expression, for Cauchy media, of the first Lame coefficient 
Xcauchy, and v+Vc corresponds to the classical second Lame coefficient pcauchy Comparing then 
equation p6| ) with equation p2| ), we find : 

Ay 2 +Az 2 n 
12 



Ax +Az 
12 
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For a given h = Ax = Ay = Az, we see that the modified equations for the scheme are those of a 
Cosserat generalized continuum, with second-order accuracy, and the coefficients verify a + (3 = 
and 7 = 2 (i+u) ^ '■ m me case °^ an anisotropic mesh size (Ax ^ Ay ^ Az), we cannot identify 
the coefficients with the isotropic Cosserat equations, due to the presence of the Laplacian operator. 
We can however find an anisotropic Cosserat model with weighted second derivatives instead of the 
Laplacian. 

One of the main characteristics of a Cosserat generalized continuum is to exhibit a characteristic 
length for the material, l c , which describes the length of the nonlocal interactions. l c is defined as : 

12 _ 7 



In our case, we see that : 

L = ^—h 
2 

l c is of the same order as the size of the particles. In an homogenization analysis framework, S. Forest, 
F. Pradel and K. Sab have shown [ 1 1 1 that when the macroscopic length of the system is fixed and the 
characteristic length l c of the Cosserat continuum tends to 0, the macroscopic behavior of the material 
is that of a Cauchy continuum. We therefore converge to a Cauchy continuum as h tends to 0. 

As a consequence, displacement £, acceleration £, rotation 9_ and acceleration of rotation 9 in 



equations (21 1 and (22 1 converge to finite macroscopic quantities. Therefore, using the equations on 
rotation, we find : 

0= ^cur!i + 0(fr 2 ) (27) 

which is the classical definition of the local rotation of a Cauchy material at order 2. Using this relation 
in the equations of displacement, we find the equations of linear elasticity for a Cauchy continuum 
medium up to error terms of order 0(h 2 ) : 

* = 2(if^ + (i + .Ki-2^ div ^ + ° {h2) 

and taking jj curl of this equation, we find the equivalent equation on rotation up to error terms of order 

0{h 2 ): I 

^^—M + 0^) (28) 



We recover a second-order accuracy on the rotation 9. As equation (27 i shows, 9 is a derivate of £, 
and we should expect only first-order accuracy using a second-order accurate method on £. We have 
therefore improved the accuracy on 9 using the Discrete Element method. 



4. Preservation of the Hamiltonian structure by the time integration scheme 
4.1. Description of the scheme 

The model built has a Hamiltonian structure. To preserve this property after time discretization, we 



use a symplectic time integration scheme. As the system (12i-(15i is a constrained Hamiltonian 



system fT5] Sec VII.5], it is natural to use the following RATTLE scheme [ 1 1 with time-step At : 
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T? +1/2 = T?-^|^(X",g") (29) 

= ^ _ At^L^n gn) + ^ngn (3()) 



= X n j + —T n T +1/2 (31) 

g„+i = gn + /± t ^n+i/2^-i (32) 

where A" is such that Q" +lT • <3™ +1 = Id (33) 

^+1 = z n+l/2 _ ^I^L^n+l gn+l) (34) 

2 uX_i — 

g+i = g+i/a _ At^y^+^gH-ij + ^a; +1 £™+\ (35) 

where is such that Q n + lT ■ PJ> +1 ■ iZ; 1 + iZ; 1 ■ I^ +lT • = Q (36) 
where A™ and A" are symmetric matrices, the Lagrange multipliers associated with the constraints 



( 33 1 and ( 36 1. We denote the scheme (29 H 36 1 by : 



{x n+i ,Q n+i ,T n+i ,E n+i ) = fA t (r,o n ,r,f) 

The proof for RATTLE'S symplecticity can be found in [30) . As a consequence, in the absence of 
exterior forces, the energy of the system is an invariant of the system, and is preserved by the numerical 
integration in time. More precisely, the error is of order 0(e~st ) over a time period of e^t , with k > 
independent from At (T3J. This yields the stability of the simulation over long time periods if the time 



step is chosen sufficiently small. In addition, we directly derive from (29l-(36i that the linear and 
angular momentum are exactly preserved. 

Another important property of the RATTLE scheme is its reversibility. Starting with the knowledge 
of positions and velocities at time (n + 1) At, we recover the positions and velocities at time nAt with 
the following scheme : 

(Q Tn ,g Rn ,P Tin ,R Rn ) = *- M (Q T<n+v Q Rtn+i ,P T ,n + i,E Rin+1 ) 

As a reversible scheme, RATTLE is of even order, and as it is consistent, it is a second-order scheme. 

RATTLE has the advantage of enforcing explicitly matrix Q n to be a rotation matrix, and at the 
same time be explicit in time. However, the nonlinearity of the constraint on Q n needs to be solved 



with an iterative algorithm, which will be addressed in section 4.3 
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4.2. Implementation with forces and torques 

For effective implementation of the RATTLE scheme, a difficulty arises from the fact that we do not 
necessarily have a direct access to ^(X", Q n ) and ^ L (X n , Q™)> as we compute the expression 

of forces and torques rather than the functional Uh- In the particular case studied here, we could 
impose directly Uh in the computation of velocity and position, but in that case, we would not be 
able to treat non-conservative exterior forces and torques, and the extension of the method to more 
complex behavior laws for the material would become unfeasible. To that end, we have chosen to 
recover |^(A n ,<3 n ) and §Q L Q£ n >Q n ) from the expression of forces and torques. We prove, in 

Appendix 



IV 



that the equations to be solved have the same form as (|29 

dU h 



35|), replacing with 



= - 2~2 JeVl Ejj and W with -kliMDQ", where M" = £ JeV , K u , and changing the 
Lagrange multipliers. 

In order to implement the scheme, without having to compute matrices A" and A , we follow once 
more (B] Sec VII.5]. We set : 

Li = Q 7 j T ■ 

gn+l/2 _ qjiT _ pn+1/2 _ 

We use the following algorithm : 

• We start the time step knowing A™, Q n , Zj 1 ' 1 / 2 and T™" 1 ^ 2 (in the first step, these last two 
elements are the null matrix and the null vector). 

• We compute the forces and torques in a submodule of the code, using only positions A" and 
Q". 

• The displacement scheme is written : 

= yn-1/2 + M ^n 
X n+1 — A™ + ^ rpn+l/2 

1 1 mj ~ 1 

• Then, we use the rotation scheme : 



- compute = Rj ■ ir 1/2 - ir 1/2 • Ut + Atc if ■ tfMj) ■ 9. 

- Find Z n T +1/2 such that : 



M + At^ +1/2 is orthogonal 

r?n+l/2 jj jj _ ^n+l/2 T _ j^n 



(37) 



- Compute Q™ +1 =Q n r -(Id + Atg> +1/2 ) 



—i —i 



We can observe that all those steps are explicit, and that the only step that requires an iterative 
resolution is (37 1. Following |T5|, we use the quaternion iterative method to solve (37 1 for Z_ n+1 ^ 2 - We 



describe that 



method in the next subsection. 
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4.3. Resolution of the nonlinear step 

Note that A™ is a skew-symmetric matrix, which can be written as : 




Equation ( [37] ) now reads : 



(38) 



To impose the second line of ( 38 I, we write the matrix Id + AtZ_™ +1 / 2 with the quaternion notation : 



M + At^ +1/2 = (e 2 , + e\ + e\ + e\)]A + 2e J| + 2E 2 



with : 



E = 






-e 3 


e2 







-ei 


-eg 


ei 






We make use of the property that every orthogonal matrix can be written in this form, and that condition 



+ eg = 1 ensures that such a matrix is orthogonal. Equation ( 37 1 is hence equivalent to 



solving for eo, ei, e 2 , e3 the following quadratic system of equations 

2{d 2 + d 3 )e ei + 2(d 2 - d 3 )e 2 e 3 
2(di + d 3 )e e 2 + 2(d 3 - di)eie 3 
2(di + d 2 )e e 3 + 2(<ii - d 2 )eie 2 



en + ef 



el 



Atai 
Ata 2 
Ata 3 
1 



(39) 



Existence and uniqueness do not hold for this set of equations. In the simple case where a± = 
ct 2 = a 3 = 0, there are distinct solutions for (eo, e\, e 2 , e 3 ) : (1, 0, 0, 0) (in that case, Z_ n+ ? = Id) , 

(0, 1, 0, 0) (in that case, Z" +5 represents the axial symmetry around axis x), (0, 0, 1, 0) (associated 
with the axial symmetry around axis y), (0, 0, 0, 1) (associated with the axial symmetry around axis z), 
and their opposites which represent the same transformation. There is a deep physical reason for that 
non-uniqueness : dynamically speaking, the rigid body is totally represented by its equivalent inertia 
ellipsoid (the ellipsoid with the same axes of inertia and moments of inertia), which is invariant under 
the axial symmetries around the inertial axes x, y and z. As the rotation Id+ AtZ^ +1 ^ 2 is an increment 
of the global rotation of the particle, we select a solution "close" to identity, in a certain sense. 

The existence and uniqueness in a neighbourhood of identity can be obtained from the equivalent 
formulation of RATTLE using the discrete Moser-Veselov scheme, with a fixed point theorem applied 
on equation (17) of reference 1 14]. We have found an explicit bound on the time-step A< for the iterative 
scheme to converge, and ensure existence and uniqueness in a neighbourhood of identity. It is derived 
in Appendix [V] We use the following iterative scheme 1 15 1 : 



We start with (eg, e", e 



1) e 2: e 3> 



(1, 0, 0, 0) (which represents identity). 
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At each iteration, we compute 



Je+l 



Ata x - 2(d 2 - d 3 )e%e% 
2(d 2 +d 3 )eg 

Ato 2 -2(d 3 -rfi)eie| 
2(di + d 3 )eg 

Ata 3 - 2(di - d 2 )ef ef 
2(di+d 2 )eg 



1 - (ef +1 )2 - (e^ 1 ) 



fc+i\ 2 



(e 3 fc+1 ) 2 



Let us introduce 



B C^-) = { ( e o,e 1 ,e 2 ,e 3 )/e^ + e 2 1 +e 2 + e 2 3 = 1, e 2 + e 2 2 + e 2 < ~ 



When the time-step At satisfies the condition : 

At 



l«i| 
h 



h 



l«3 

h 



< 



0.26 



(40) 
(41) 
(42) 
(43) 



(44) 



the algorithm (40 1-(43 I converges with a geometrical speed to the unique solution in B(^). 

Let us observe that li and D scale as ph 5 . In addition, as P_ l = j{^lj)Q D_ T , is of the order 

of HO/11- Using the expressions (18i and (19 1, and the fact that a n , a s and at scale as h 2 , we obtain 
that Mj is of the order of Eh 3 . Condition (44 > therefore gives us a constraint on the time-step of the 
following type : 

At 2 E 



Ai ft, 



P 



< C 



(45) 



where C is a constant. This is the natural CFL condition for an explicit scheme on rotation, with 
the typical celerity of the compression and shear waves in the material. 



5. Numerical results 

In this section, we present several challenging test cases. First, we address Lamb's problem, which 
allows us to examine numerically the precision of the method in the case of small displacements against 
a semi-analytic solution. The presence of surface waves is the most difficult part of the problem, and 
the results appear to be satisfactory. We examine the conservation of energy on the case of a three- 
dimensional cylinder submitted to large displacement. In the end, we also demonstrate the ability of 
the method to tackle static rod and shell problems using the same formulation, on the cases of the 
bending of a rod and of the loading of a hemispherical shell. 

5.1. Lamb's problem 



We have simulated Lamb's problem (see \21\) : a semi-infinite plane is described by a rectangular 



domain, with a free surface on the upper side, and absorbing conditions on the other sides. On a surface 
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particle, we apply a vertical force, whose time evolution is described by a Ricker function (the second 
derivative of a Gaussian function). We observe the propagation of three waves : inside the domain, a 
compression wave of type P and a shear wave of type S, and on the surface, a Rayleigh wave. We also 
have a P-S wave linking the P and the S waves, which is a conversion of the P wave into an S wave after 
reflection at the surface. In the case of a two-dimensional problem, the intensity of P and S waves is 
inversely proportional to the distance to the source, and the intensity of the Rayleigh wave is preserved 
throughout its propagation. 

We have chosen the following characteristics for the material : the density is p = 2200 kg.m~ 3 , the 
Poisson coefficient is v = 0.25, Young's modulus is E = 1.88.10 10 Pa. The velocity of P waves is 
therefore approximately 3202 m.s -1 and the velocity of S waves is 1849 m.s -1 . 

The force applied is a Ricker of central frequency 14.5 Hz, that is, with maximal frequency around 
40 Hz. The minimal wave length for P waves is therefore 80 m, and the minimal wave length for S 
waves is approximately 50 m. In the rest of this subsection, we call "wave length" this minimal wave 
length of 50 m. We indicate the discretization step in terms of number of elements per wave length. 

Lamb's problem has the interesting particularity of having a semi-analytic solution : Cagniard's 
method is described in |6j. We have compared our results with this exact solution and thus estimate the 
numerical error of the scheme. The comparison between the numerical results and the semi-analytic 
solution obtained at 300 meters from the source, on the surface, with Ax — Ay = 5 m (10 points per 
wave length), is shown on figure [2] 




(a) Horizontal displacement 




(b) Vertical displacement 



Figure 2. Displacement at the surface, 300 meters from source, with Ax = 5 m, Ay = 5m (10 points per wave 

length) 

We compute the same result with different spatial discretizations, with Ax = Ay. As expected, 
refining the spatial discretization decreases the error. The velocity of the different waves agrees with 
the exact solution, and the amplitude of the waves is accurately captured with more than 10 elements 
per wavelength. The accuracy of the method cannot compare with that of spectral elements (5 points 
per wave length), but it gives better results than classic second-order finite elements (30 points per 
wave length), and mostly on the surface, where we recover the non-dissipative Rayleigh wave. This 
is probably due to the introduction of parameter 9 which helps us simulate the rotation of the particle 
precisely, instead of recovering it as a Taylor development of the displacement, thus losing one order 
of accuracy for rotation. 

If we measure the L°°-error on vertical displacement at 300 meters from the source, with an angle 
of 60° with the horizontal axis, we obtain an approximate slope of 2 fitting the points (figure B). This 



confirms the results of subsection 3.1 as to the second-order nature of the spatial scheme. 
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Discretization step (m) 

Figure 3. Linear fitting of the log-log diagram for the numerical error against the spatial discretization step 
5.2. Conservation of energy 

In order to illustrate the conservation of energy by the scheme, we model the evolution of a pinched 
cylinder. The cylinder has a radius of lm, a height of 2m and a width of 1cm. The physical 
characteristics are that of steel (E = 210000 MPa, v = 0.25). The cylinder is discretized with 50 
elements on the perimeter, 20 elements on the height and one element in width. Opposite forces are 
applied on two sides of the cylinder, pinching it. At the initial time, the forces are removed, and the 
cylinder is left free. We simulate the system over 500,000 time-steps, corresponding to 45 oscillations 
of the first mode of the cylinder. The large number of time-steps required reflects the fact that a number 
of smaller local oscillations propagate at high velocities, and that the cylinder is very thin. On figure|4] 
we observe an excellent preservation of the energy. The configuration of the cylinder at the moment of 
release is shown on figure [5] 

The preservation of energy is quite satisfactory, even with large displacements in a three-dimensional 
geometry. 



5.3. Static shell test cases 

In order to show the versatility of the method, we compare the static deformation obtained with Mka3D 
(adding damping to the model) to the second and fourth benchmarks for geometric nonlinear shells 
found in pTj . 

The first benchmark considered is that of the cantilever subjected to an end moment A4. Let N be the 
number of discrete elements in the length of the beam. We take one element in the two other directions. 
We immediately see that at the equilibrium, for each particle /, the sum of forces is null, and using the 
boundary conditions, the force Fjj between particles is always null. The sum of moments is also zero, 
and is equal to the end moment A4. As Fjj = 0, if we denote 6n the angle between two consecutive 



particles, using (51 



M u = M{j =^rsm6 N (46) 



If we take the maximum end moment M. m ax — 2ir^-, which is the theoretical moment applied to 
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Figure 5. Initial configuration of the cylinder 



bend the beam into a circle, we obtain : 

N9 N = N arcsin ^ J (47) 

As N tends to infinity, the deflection angle of the end N6n tends to 2ir with second order precision, 
which indicates a second order convergence to the theoretical solution. This convergence has been 
checked in practice. 
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The second benchmark considered is a hemispherical shell with an 18° circular cutout at its pole, 
loaded by alternating radial point forces T at 90° intervals. The shell is discretized by 16 elements in 
latitude, 64 elements in longitude and one element in thickness. The initial and deformed geometries 
are shown on figure [6] The radial deflections at the points of loading A and B are compared with the 
results obtained in pi) in figure [7] Our results are in very good agreement with the benchmark. 




Figure 6. Initial geometry and deformed geometry at T = 4007V for the hemispherical shell subjected to 

alternating radial forces 




Figure 7. Load-deflection curves for the hemisphere shell at the loading points A (left) and B (right) 
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6. Conclusion 



In this paper, we proposed a numerical discretization of material continuum, allowing for the simulation 
of three-dimensional wave propagation as well as shell or multibody dynamics, in a monolithic 
way. It is consistent with the equations of elastodynamics at order 2 in space and in time, and we 
numerically recover the propagation of seismic waves in the body of the material and at the free 
surface. Furthermore, the dynamics of the system are written in the form of a Hamiltonian dynamics. 
Using symplectic schemes, we correctly reproduce the preservation of the system energy. This ensures 
numerical L 2 -stability of the scheme, and allows long-time stable simulations with large displacements 
and large deformations. As the method is entirely local and requires no matrix inversion, it can 
be easily parallelized with domain decomposition. The main restriction is the size of the time-step 
due to the explicit nature of the integration scheme. This could be remedied by using asynchronous 
symplectic integrators in order to have local time refinement at small elements and a global larger 
time-step [31]. This work can be seen as a first step towards using more complex constitutive laws 
(while still maintaining stability of the scheme), and towards coupling particle dynamics simulation 
with a fluid dynamics simulation for fluid-structure interaction. 



ACKNOWLEDGEMENTS 



The first author acknowledges the support of CEA under Grant n° 1045. 

We would like to thank Serge Piperno, Tony Lelievre, Frederic Legoll and Eric Cances (Cermics and UR Navier, 
Ecole des Ponts) for useful discussions and advice on the mathematical and computational aspects of this paper. 
We also thank Karam Sab (UR Navier, Ecole des Ponts) for pointing us the similarity of our model with Cosserat 
models. Thanks are also due to Gilles Vilmart for discussions on the resolution of the quaternion scheme. 



APPENDIX 

I. Expression of the equivalent volumetric deformation with a free surface 

We need to account for the boundary condition rr • n = at every free surface of the particles. We have seen 
that the discrete equivalent for g_ ■ n is Fj j. For a given particle /, we assume that the particle is 



2.3 



in section _ _ _ 

surrounded by real particles J G Vi, and by 'ghost' particles J S V\ at every free boundary. The position of these 
particles is ajusted in order to satisfy the boundary condition. 

The equivalent deformation of particle / can be expressed as in the bulk of the material : 

„ ^ 1 SlJ . TT-^ 1 Su . 

£l = 2^ 2 "U~— IJ ' + 2 ~V"— IJ ' - IJ 
Jevj 1 Jev } 1 

For a ghost particle J € V}, the boundary condition F_j 7 ■ n t 7 = boils down to : 

W^~ IJ ' - IJ + g " + = (48) 



Summing 1 48 1 over the ghost particles, and using the fact that the free volume V\ satisfies 

Jev 1 , 
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we find that the deformation of the links with the ghost particles should follow the equation : 

v- lSu_ 3v Vj lgrj 

2- 2 Vj — IJ ' ~ IJ - 1 - 21/ Vj + j%,V} fa 2 V/ — " ^ - /J 

Inserting this relation in the expression of e", we check that : 



Jev 7 v/ ^ 1-21/ 



We denote : 



II. Expression of the coefficients for the flexion and torsion of the particle links 



I'u= II (XPu ■ SjjfdX 

J J S T J 

I\j= If {XPu-tjjfdX 

J J St j 



(49) 
(50) 



the principal moments of the interface between particles I and J, we require that : 



a„ + a s 



a„ + at — 



a s + a t 



Elfj 

Elfj 
Su 

Ejlfj + Ijj) 
2(1 + v)Su 



(51) 



The expression of the a is given by : 



_ (l + 2v)E t , 

Un ~ 4(1 + v)sJ Iij+Iij) 



a s = 



E 



4{l + u)Su 
E 



((3 + 2v)Iij - (1 + 2iy)ljj) 



at = ((3 + 2v)IlJ - (1 + 2v)I ' j) 

4(1 + v)bu 



(52) 
(53) 
(54) 
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III. Derivation of the forces and torques from the potential energies 

The derivation of potential energies is straightforward : 

dU t ^ Su E 



dXj 



.rev, 

dUg _ 

- 1 Jev, 



®Ut \ ^ 57./ E (, 



^ fi + „Ki-2^ 5lje?J ^ J ® 



Sg, ^(1 + ^(1-2*,) 

at/ 

<9Q 



=/ Jevj 

Using the expression of the force F u between particles I and J : 

Su E 



— IJ ~ D° 



E . „ £V „ / 1 A 1 , A \ 

TT^" + 5/J (i + «,)(i- 2,) £ " + b77^ j " i>7; ( — " ' - Ij) - IJ ) 



we obtain : 

m i*i =ti = F u 
For the rotational part, it is easily obtained that : 

m^i) = m^R - qmi) = ■ £ T - % ■ e! 

Deriving in time, we obtain : 



^»-M-(f)a/ + a,(f)' 



Using the fact that : 



(a <g> 6) • Q = a <g> (Q 1 • 6) 

we get : 

dQ 



m h • £/ = - £ §^ if^*» ® (2, ' ^) (55) 



^ • 2/ = - £ ( 1 + ,) (1 - 2|/ ) ^W « (Q • xfe) (56) 

+at(Q j -t I j)®(Q zI -tij)) (57) 
Denoting symm() and skewQ the symmetric and skew-symmetric parts of a matrix, we note that for any a and 

b: 

j(a A b) — — skew(a (g) 6) 
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Using the expression of the torques M \ ; and M_{ 7 : 

= §r ~ (£, • m±) a am,, + (1 + ^_ 2t/) ^(g • m±) a 

Mjj = -§r(a»(Q -n?j)A(Q -n°j) +a s (Q ■ s u ) A (Q +a t (Q ■ t u ) A (Q 

equation l |15^ gives us the equation on the angular velocity : 

d 
dt 



R^Qi) = E mL + mL 

JGV 7 



IV. Details on the implementation of the RATTLE scheme with forces and torques 
For forces, the relation is simple : 

H<*,g) = -££,, 

For torques, we have : 

FIT I 

where is the symmetric matrix of Lagrange multipliers associated with constraint Q ■ Q 1 — Id. On the other 
hand, 

l ( E Hu J -£, ■ 3/ + £,-£,'- 2, - £, T - 2, ' £, T 

as the A^ are symmetric. Therefore, there exists a symmetric matrix A J such that : 

dU h , 



o^9) =[-2 



\jevr ) , 



We denote : 



& = E Ejj 

.lev, 

m = E m u 

JGV, 

where forces F_ } 7 and torques M_j 7 have been computed with positions X™ and Q n . 
We can rewrite equations \29\ to {35} as follows : 

yn+l/2 = £ „ + At^n (5g) 
^n+1/2 = + ^|(M")£" + ^ (4l + 4;' )£" (59) 
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gn+1 = x« + ^r" +1/2 (60) 

m; 

gn + 1 = gn + ^ 

where A_" is such that Q" +lT • Q" +1 = M (62) 

yn + l = yn + 1/2 + ^ £ n + l (63) 

= £; +i/2 + ^j(m? +1 )£; +1 + + a; +1,0 )£; +1 , (64) 

where A n is such that Q n+lT ■ P n+1 ■ D~ l + C" 1 ■ P n+lT • Q" +1 = (65) 
— i =i —i — i — i —i =i - 



V. Resolution of the nonlinear step of the RATTLE time-scheme 

In this appendix, we examine the resolution of the nonlinear step of the RATTLE time-scheme described in section 
|4.3| We determine conditions on the time-step At that ensure convergence of the iterative algorithm <|40[> — (|43[> 
in a certain neighbourhood of identity, and we conclude on the existence and uniqueness of a solution in this 
neighbourhood. 

We denote 6(0, r) the ball of center and radius r : 

B{0,r) = {(ei,e 2 ,e 3 )/e 2 + e 2 . + e 2 < r 2 } 

Using the numerical scheme described in section 4.3 we first show that it stabilizes a ball included in B(0, 
under a CFL-type condition on At. We then show convergence in that same ball, and we conclude on convergence 
to the unique fixed point. 

V.l. The iterative scheme is bounded 

Starting with a given (go, ei, e%, e 3 ) computed in the previous iteration, such that e 2 , + e 2 + e| + e| = 1, the 
iterative scheme l|40|-i 43 i gives the new quadruplet [e% , e\ , e% , e 3 ) defined by : 



ei 



e 2 



Ata\ —2{d,2 — d- i )e2e- j , 

2(d 2 + d 3 )e 
Ato 2 - 2(rf 3 - di)eie 3 
2(di + d 3 )e 
« _ Ata3 — 2(di — d2)eie2 
63 ~~ 2(di + d 2 )e 

el = y/1 ~ (el) 2 ~ (4) 2 - (e 3 *) 2 

For this scheme to be well-defined, (e*, e|, e 3 ) should be in 6(0, 1). We impose a stronger condition, with 
(ei, 62, e 3 ) and (e£, e^e^) in ,6(0, /3) where /3 is less than |. 
Suppose that : 

e 2 + e 2 + e 2 < 

We want to have : 

(elf + (e 2 f + (e 3 ) 2 < P 
As e 2 , + ef + e| + e 2 = 1, we also have e 2 > 1 — p. Since : 
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we obtain : 



|ei| < 



Let us define I\ = d% + ^3, I2 = di + c?3, ^3 = di + 0(2 and 



(| Attn I +/9|d 2 -d 3 |) 



/(/?) 



1 



4(1-/3) 



A< 2 / jojf ]o&£ ]o«f 



-2/3 At 



n 



T 2 



1 .i 



+P 1 



n 



+ 



11 



+ 



11 



then the previous assumptions imply that : 

(el) 2 + (e^) 2 + {e 3 f < /(/?) 
Therefore, a sufficient condition for the scheme to be bounded is /(/3) < /3. We know that : 

I d 2 — cfe I M2 — dz I 



/1 



o?2 + dz 



< 1 



as the di are positive. Then 

™ * 4(1^3) 



Ar 



Mil 2 , M2I 2 , Msl 2 



n 



T2 
1 > 



I 2 



+2/3At 



Qi , |Q2| |« 3 | 

/l 1 2 / 3 



+ 3/3' 



Hence, a sufficient condition for /(/3) < /3 to hold is : 



At 



1? i! 



Let us define : 



1 { !'H| 2 M2I 2 , M3| 2 
+ r o + T 2 
J 3 



B 

C 



+ 2/3 At 



\oi2\ |Q3 
T T T 



+ 7/3 - 4/3 < 



|Ql |a 2 | |<*3| 
111 

I 1 2 1 1 2 1 1 2 

|ai| , |«2| |a 3 | 



if 



r2 

J 2 



I 2 
J 3 



A sufficient condition to obtain 1 66 1 is to have At < At with : 

-2/35 + ^/4f3 2 B 2 - 4(7/3 2 - 4/3) C 



At = 



2C 



(66) 



As we supposed that < /3 < | < =, 7/3 - 4/3 < 0. We also know that B < 3C and C < B 2 , and it follows 
that : 



ft > 



2^/ ^-/3 



In the end, we have the following lemma : 
Lemma V.l. Let us choose < /3 < | anrf At > jmc/; /ftaf : 



At 



Ql 1 "2 1 Q 3 

III 



< 2 



(67) 



//(ei,e 2 ,e 3 ) € 23(0, S),then {e\, e* 2 , e* 3 ) £ B(0, V/3). 
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V.2. The iterative scheme is a contraction 

Following the previous subsection, suppose that (ei,e2,e3) and {fi, fz, fs) are in £?(0, \//3), and let eo = 
\J\ — e\ — e| — e§ and fo = \/l — f( — /f — /| . We define e* and /* as before. We show here that 
||e* — /*|| < p||e — /||. withO < p < 1. 
We compute : 



f * (d 2 ~ ^3) 



he 

2- 
ii 



(/2 - e 2 ) — - — +(/ 3 - e 3 ) 



. fo — eo .» 

H — : — ■/ 1 

eo 



We then use the fact that < 1. As the same type of results hold with a circular permutation of indices 



x, y and 2, we let ||- 1 the euclidian norm in K on (ei , e2, 63), and we find : 

lie* - /* f < 2 ^ 2 j ±A 2 j (/l - ei ) 2 + 2 l 2 j y 2 j (/ 2 - e 2 ) 2 

+ 2 ( ^ )2 e + 2 ( ^ )2 (/3 - 63)2 • > - es) (^) 

4,. w. > / fi + ei \ / fj + eg \ 4 , f fi + ei\ f f 2 + e 2 



2 7 V 2 7 e 2 w /w ' V 2 



Since : 



+ 2 (/D 2 + (/y + (/3) 2 (/o _ eo)2 



\ 2 / \ 2 

/f „ ^ / 72 + e 2 \ , . s2 / 73 + e 3 
(/2 - e 2 ) +(/ 3 - e 3 ) 



we have : 

lie* - /* || 2 < 1 (ll^flle - /|| 2 + |iri| 2 (/o - e„) 2 

We also have : 



(/o-e ) 2 <i^||e-/|| 2 



In the end, we obtain the upper bound : 



e+f ||2 / 11 -r* 11 2 



If we take the same hypotheses as in the first subsection, that is, (ei, e 2 , 63) € B(0, v7?) and (fi, / 2 , fi) £ 
B(0, v^), and ft such that (et,e* 2 ,e* 3 ) € B(0, V?) and (/i*,/ 2 *,/ 3 ) € B(0, v^), then due to the convexity of 
B(0, y^, we have : 

ll^f</3 



and moreover, as eg > 1 - /? et /„ > 1 - f3, then (^4^) > 1 



Then : 



2 



eo+/o ) 2 J ~ 1-/3V 1-/9/ (1-/?) 



In order to have a scheme which is a contraction, it is sufficient to impose : 

^ <1 
(1-/3) 2 " 

As < fj < i, it is sufficient to choose : 

8 < 2- \/3 
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V.3. Optimization on constant f3 

Optimizing the stability condition \(fl) on At, we obtain the following optimal value of /3 : 

= L=^E w . 17 

V.4. Conclusion 

If we take the time-step At such that : 

M (M + M + M\< j p max -^ _ fimax « o.26 

\ il J2 i3 / V o 

then the iterative scheme starting with (1, 0, 0, 0) converges to the unique solution of the nonlinear problem in 



B(0, J 7 i4^ )» an d the convergence speed is geometric with a rate p < 1. In addition, p < 28 — 6\/2T w 0.5. 
We thus have proved existence and uniqueness of the solution in B(0, ^). 
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